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Motivated by the decade-long debate over the issue of criticality supposedly observed 
in nuclear multifragmentation, we propose a dynamical lattice model to simulate the 
phenomenon. Its Ising Hamiltonian mimics a short range attractive interaction which 
competes with a thermal-like dissipative process. The results here presented, generated 
through an event-by-event analysis, are in agreement with both experiment and those 
produced by a percolative (non dynamical) model. 

1. Introduction 



In the early 1980's, the experimental observation of power law-like distributio 
of fragment yields in nuclear reactions induced by high energy proton projectiles u\ 
initiated an ongoing controversy on the origins of this behavior. When a proton or a 
heavy ion collides with a heavy nucleus, the target can break into a number of nuclei 
with atomic number Z > 2 in addition to several nucleons and alpha particles. The 
number of heavy fragments induces a classification of these processes into three 
different regimes: fission, with two heavy fragments; spallation, yielding one heavy 
fragment and a few light ones; and multifragmentation, where the product includes 
no heavy fragments, but is composed of several light nuclei of varying sizes 0. There 
is general agreement that the boundary between light and heavy fragments lies 
around Z = 20 □. The processes of fission and spallation are quite well understood, 
but such is not the same with multifragmentation. One of the reasons for this 
lack of knowledge stems from experimental difficulties - a detailed investigation of 
fragmentation requires coincident measurements of the multiple fragments formed. 
Most of the experiments reported in the literature measure little more than the 
inclusive mass yield of fragments - and so cannot make a clear distinction between 
those that correspond to each kind of phenomenon involved. Most of the exclusive 
experiments, on the other hand, are emulsion experiments with inverse kinematics a, 
and suffer the drawback of poor statistics. The lack of reliable coincidence data has 
forced theoretical investigations to concentrate on the explanation of the mass yield 
curve. 
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For Z < 20, it has been pointed out that this curve is compatible with a power- 
law dependence for the cross-section a ~ Z~ T ', where r ~ 2.2 is practically indepen- 
dent of the beam energy and the exact composition of projectile and target. This 
parametrization of the mass yield curve, along with the measured value of the ex- 
ponent, is precisely what is expected for the transition between liquid and gaseous 
phases in nuclear matter □ - this transition takes its name frojji a parallel with the 
process of droplet formation in a similar transition for fluids U. This fact has been 
pointed out as a clear indication that multifragmentation is strongly related to that 
phase transition. 

The currently accepted scenario for the fragmentation of nuclei begins with a 
quick isentropic expansion, as the target nucleus is hit Q. The resulting compound 
nucleus has its density diminished in this process, and the fragmentation begins 
with a cracking of the system. In this way, primary fragments are formed, still 
highly excited. Deexcitation of these fragments is acompanied by the formation of 
secondary fragments 13. During the whole process, the system has enough time for 
thermal equilibration; we have two different time-scales involved - namely, in this 
case, we have fast formation of cracks and slow thermalization -, one of the known 
conditions for the installation of a complex regime. The theoretical understanding 
of this process is very difficult, since we are in the presence of a true many-body 
problem, due to long range correlations between the clusters that are being formed. 

Among the theoretical models that have been proposed to study this phe- 
nomenon ( for a brief review, see Ref.p] ), by far the most popular one is a bond- 
percolation model on a cubic lattice @. It has been shown that, independently of the 
knowledge of the parameter that would control the approach to criticality, a number 
of successful comparisons could be made between predictions of this model and ex- 
perimental results £3. By studying the moments of the cluster size distributions on 

a event-by-event analysis, and assuming a scaling property for these distributions 
, one can show that a strong correlation between these moments must be present 
in the vicinity of the critical region. Nevertheless, the purely geometrical and static 
character of such a model brings about the need for one that could rely on some 
dynamical interaction between the nucleons participating in the process. 

In the next section, we present a dynamical dissipative lattice model that re- 
sponds to that need. Some details of its computer implementation are then dis- 
cussed. The last section is dedicated to a presentation of our results, and to a 
comparison between them and those coming from the percolative model. 

2. The Model 

Initial versions of the model we present in this paper were used in the context of 
surface wetting E3 and of drop formation on the leaky faucet problem liSo. Its use 
in nuclear collisions was foreseen in the analysis of residual mass in the evaporation 
of hot nuclei Ej. The success obtained in that application encouraged us to the 
extension presented here. 
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The nucleus, in its initial configuration, is a dense cluster of occupied cells in an 
otherwise empty cubic lattice - using the terminology of magnetic systems in the 
context of the gas-lattice model, a cluster of spins up surrounded by spins down - 
subject to an Ising-type interaction given by the model Hamiltonian 

n = -j SiSj, (i) 

<i,3> 

where the summation extends to nearest and ncxt-ncarcst neighbours and Sij = ±1. 
With this last prescription, we intend to mimic a surface tension, having in mind a 
liquid drop model for the nucleus. The initial excitation energy of this compound 
nucleus - target plus absorbed projectile - is associated to a temperature parameter, 
which will control the transition probability between two different configurations of 
the system. The nucleus is then subjected to a dissipative Metropolis dynamics 0, 
generating a Markovian chain of configurations. This dynamics involves: 

(i) a random choice of an occupied site at the cluster perimeter, followed by a ran- 
dom choice of an unoccupied site, not contiguous to the first, also at the perimeter; 

(ii) a double flip of the spins at these two sites, decided by the Metropolis rule - 
thus involving the temperature parameter already mentioned; 

(iii) if the double flip occurs, it is verified if this configuration is still a connected set 
of spins up; if it is not, we have the formation of a primary fragment, whose mass 
is contabilized, and energy is dissipated, in the form of a decrease in temperature 
(intermittent, fast cracking energy dissipation); 

(iv) in either case, an additional small decrement in temperature is promoted, sim- 
ulating the emission of radiation by the system (continuous, slow thermalization). 

The recently formed fragment is erased, and this dynamics continues until a low 
temperature is attained, when no more fragments would be formed. 

The distribution of fragments thus obtained is characterized by its moments. 
If rii(s) is the number of fragments with mass s obtained in the event i, its k th 
moment is computed by 

Mi(k) =^ S fe n, ; ( S ). 

s 

In a normal thermodynamic system, these moments would diverge in the thermo- 
dynamic limit at the critical point, for k > 1, with critical exponents 

M(k,e) ~ e-» k , 

where e measurea^thfi distance to the critical point. If the distribution n(s) has the 
scaling property B'tll 

n(s,e) ~s- T f(es°), 

where r and a are two critical exponents, one can show that moments of differ- 
ent orders must be correlated. The relation betw een the exponents of moments 
divergence and those of the distribution is given by t£l 

fi k = -(r- 1 - k)/a. 
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Since we are working with intrinsically finite systems - the nuclei - the moments 
will remain finite, even for k > 1. The normal signature of critical behavior - the 
divergence of the moments of the distribution - is washed out by the finiteness of the 
system. This is one of the difficulties that faces theoretical work in small systems. 
Nevertheless one can use the correlation between moments of different orders to 
examine the surviving traces of critical behavior. 

It is more natural to work with normalised moments 

Si{k) = Mi(A)/Afi(l) 

instead of the regular ones. Then it can be shown that 

log(5i(3)) = A3 /2 log(5 2 (2)) 

and 

log(5 t -(5)) = A 5/2 log(5 l (2)) 
where, with the usual identification 7 = [12, 

A3/2 = 1 + l/i77. 

and 

A5/2 = 1 + 3/(77. 

3. Computer Strategy 

We will adopt the C programming language sintax in this section, with some 
minor explanations. Our strategy follows two rules: storing data (spin states) bit 
by bit on computer words; and treating them almost exclusively by bitwise logical 
operators. Obviously, this saves computer memory by a factor of 32 (for 32 bit word 
computers). Better yet, computer time can be saved by a similar factor whenever 
parallel updating is possible; such is the case in some of the routines discussed below. 
Many general tricks designed to implement this strategy, similar to the particular 
ones used here, are discussed in detail in references [S'El]. 

Our model nucleus resides on a 32x32x32 cubic lattice. This lattice is mapped 
onto a vector L[r] (r = 0,1,2... lastword), that keeps the current nucleus shape, 
through the rule: element (i,j,k) is represented by the k th bit of L[j * 32 + i] 
(lastword is a mnemonic for 32 * 32 — 1). The rationale behind this choice lies in 
the fact that, although a mathematical operation is needed whenever a particular 
cell must be accessed, that is exactly what is done at machine code level - so very 
little extra time is spent. On the other hand, most of the really time-consuming 
operations are done on a sequential basis; in these cases, storing the array as a 
single-index vector avoids the need to go through multiple accesses by pointers, 
which is the way arrays are usually implemented on computers. 
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As an example of all that has just been said, we write down the algoritm used 
to determine the inner nucleus boundary. This boundary is stored on array I[r] by 

for(r = 0; r <= lastword; r + +) { 

I[r] = (~L[r]«l)! (~L[r]»l)', ~L[r-l] \ ~L[r+l] \ ~L[r + 32] ', ~L[r-32]; 
/* first neighbors */ 

l[r] = I[r] ; ~L[r - 1 - 32] ; ~L[r -1 + 32] ; ~L[r + 1-32] ; ~L[r +1 + 32] ; 
(~L[r-l]«l)', (~L[r-l]»l)', (~L[r+l]«l| (~L[r+l]»l) ', (~L[r-32] «1) 
; (~L[r - 32] »1) ; (~L[r + 32] «1) ; (~L[r + 32] »1); /* and second neighbors 
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I[r] =I[r] &L[r]; } 

where ~, and & represent respectively the bitwise logical operations NOT, OR and 
AND. A similar procedure gives the outer drop boundary. A pair of non-adjacent 
positions are randomly selected, one from each of these boundaries, as candidates for 
a mass-conserving, Kawasaki-like updating. The energy variation involved in this 
exchange of positions is also calculated through the use of bitwise logical operations. 

The determination of the connectedness of the nucleus after the spin updating 
is a dificult problem, due to its non-locality. There is no clear shortcut solution 
in the literature for three dimensions, and here resides the most time-consuming 
algorithm in the program. We chose the so-called burning algorithm t3, designed 
to determine connected clusters of lattice sites sharing some property. A computer 
implementation for a similar problem in two dimensions was published in Ref. |ia] . 

4. The Results 

For a close examination of the correlation between normalised moments of dif- 
ferent orders, we made log-log plots of SIS) x S(2) and 5(5) x S(2), following the 
suggestions originally made by Campi 113. In fig. 1 and 2 we show the evolution 
of these correlations with the initial excitation energy of the system. It is clear 
from these graphs that the model reproduces the strong correlation observed exper- 
imentally, as also did the percolation model. The points making an apparent arc of 
circle near the origin for low excitation energies in both instances are themselves also 
present in the experiments, and are probably related to very 'gentle' events, yielding 
mostly fragments of small mass (< 4). The measured slopes A 3 / 2 = 2.28 + 0.01 and 
A5/2 = 4.75 ± 0.04 agree with both experiment and percolation model E9. 

Observation of the data gives an indication that our model is probably more 
'selective' in the energy involved in the reaction then the experiments they are 
being compared with - the data used was from inclusive collisions El. A selection of 
the initial temperature parameter can bring the system closer to the critical region. 
This makes it a candidate for modelling the new exclusive experiments now being 
made. 

The determination of the proximity to the critical region made use of the so 
called 'Campi scatter plot' (fig. 3). This is a log-log plot of the greatest fragment 



6 J.S.Sd Martins and P.M.C.de Oliveira 



1E+2 



ft 



1E+1 — 



1E+0 



T = 10 




1E+0 



1E+3 — 



1E+2 — 



T = 40 J- 



1E+1 -= 




1E+0 



1E+2 



Q2 1E+1 1E +o 1E+1 S2 

Fig. 1. Plot of log(S'3) X log(S'2)- The system used in simulations is a 6x6x6 cubic cluster p-tptal 
number of "nucleons" is 216, permitting close comparisons with the results quoted by Rcf.[E3]. 
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Fig. 2. Plot of log(>S5) X log(5'2). The points closer to the origin on the graph on the left 
correspond to highly undercritical events, as explained in the text. 



in each event against the second normalized moment, averaged over ^vents with 
the same M(l). In this plot one can identify three different regions E3: a region 
with negative slope, corresponding to events with small energy - or undercritical 
-, a region with positive slope, or overcritical and a region where there is a great 
dispersion in the second moment. We understand that this is the critical regime we 
are looking for. In fig. 3 we show the evolution of this plot with the increase in the 
temperature parameter. The simulation with T — 40 was chosen as a representative 
of the experimental situation, and determined the temperature where the slopes A 3 / 2 
and A5/2 were measured. The measured slope of the scatter plot in the undercritical 
region, although very inaccurate, is also coherent with the result of the percolative 
model. 

In summary, we showed that a dissipative dynamic model can reproduce much 
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Fig. 3. Evolution of Campi's Scatter Plot with the initial value of the temperature parameter. As 
this value is increased, the system goes more and more into the critical regime, loosely associated 
with events represented by points at the intersection of the negative and positive slopes portions 
of the graph. 



of the observed behaviour of the distribution of fragments produced in nuclear 
reactions induced by collisions with high energy protons or heavy ions. These re- 
sults allow us to conclude that a few model dynamical ingredients, namely a short 
range nucleon-nucleon attraction, fast energy dissipation during each crack forma- 
tion and continuous, slow thermalization are enough: a detailed description of the 
microscopic interactions is not needed. A complete characterisation of the fragment 
distribution generated through the use of this model is now in progress. 
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